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Abstract 

We explore a recently proposed locally resonant granular system bearing harmonic 
internal resonators in a chain of beads interacting via Hertzian elastic contacts. In this 
system, we propose the existence of two types of configurations: (a) small-amplitude 
periodic traveling waves and (b) dark-breather solutions, i.e., exponentially localized, 
time periodic states mounted on top of a non-vanishing background. We also identify 
conditions under which the system admits long-lived bright breather solutions. Our re¬ 
sults are obtained by means of an asymptotic reduction to a suitably modified version 
of the so-called discrete p-Schrodinger (DpS) equation, which is established as con- 
trollably approximating the solutions of the original system for large but finite times 
(under suitable assumptions on the solution amplitude and the resonator mass). The 
findings are also corroborated by detailed numerical computations. A remarkable fea¬ 
ture distinguishing our results from other settings where dark breathers are observed 
is the complete absence of precompression in the system, i.e., the absence of a linear 
spectral band. 


1 Introduction 

Granular materials, tightly packed aggregates of particles that deform elastically when in 
contact with each other, provide a natural setting for the study of nonlinear waves. Under 
certain assumptions, dynamics of granular crystals is governed by Hertzian contact inter¬ 
actions of the particles m This leads to the emergence of a wide variety of nonlinear 
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waves. Among them, arguably, the most prototypical ones are traveling waves [l|-;3], shock 
waves (4j[5j and exponentially localized (in space), periodic (in time) states that are referred 
to as discrete breathers (6]-8|. 

Discrete breathers constitute a generic excitation that emerges in a wide variety of systems 
and has been thoroughly reviewed (9,101. Discrete breathers can be divided into two distinct 
types, which are often referred to as bright and dark breathers. Bright breathers have tails 
in relative displacement decaying to zero and are known to exist in dimer (or more generally 
heterogeneous) granular chains with precompression [6,11,12 , monatomic granular cl 
with defects [13] (see also (l4l) and in Hertzian chains with harmonic onsite potential 
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17 . Dark breathers are spatially modulated standing waves whose amplitude is constant at 
infinity and vanishes at the center of the chain (see Fig. [9| for example). Their existence, 
stability and bifurcation structures have been studied in a homogeneous granular chain under 
Recently, experimental investigations utilizing laser Doppler vibrometry 
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precompression 

have systematically revealed the existence of such states in damped, driven granular chains 
in |7|. However, to the best of our knowledge, dark breathers have not been identified in a 
monatomic granular chain without precompression. 

In this work, we focus on a recent, yet already emerging as particularly interesting, 
modification of the standard granular chain, namely the so-called locally resonant granular 
chain. The latter belongs to a new type of granular “metamaterial” that has additional 
degrees of freedom and exhibits a very rich nonlinear dynamic behavior. In particular, in 
these systems it is possible to engineer tunable band gaps, as well as to potentially utilize 
them for shock absorption and vibration mitigation. Such metamaterials have been recently 
designed and experimentally tested in the form of chains of spherical beads with internal 
linear resonators inside the primary beads (mass-in-mass chain) [19], granular chains with 
external ring resonators attached to the beads (mass-with-mass chain) [201 (see also |2l]) 
and woodpile phononic crystals consisting of vertically stacked slender cylindrical rods in 
orthogonal contact 22 . An intriguing feature that has already been reported in such systems 
is the presence of weakly nonlinear solitary waves or nanoptera [23] (see also [24] for more 


detailed numerical results). Under certain conditions, each of these systems can be described 
by a granular chain with a secondary mass attached to each bead in the chain by a linear 
spring. The attached linear oscillator has the natural frequency of the internal resonator in 
the mass-in-mass chain (Fig. 1 in [l9j), the piston normal vibration mode of the ring resonator 
attached to each bead in the mass-with-mass system (Fig. 9 in [201) or the primary bending 
vibration mode of the cylindrical rods in the woodpile setup (Fig. 1 in [221). 

One of the particularly appealing characteristics of a locally resonant granular chain of 
this type is the fact that it possesses a number of special case limits that have previously been 
studied. More specifically, in the limit when the ratio of secondary to primary masses tends 
to zero, our model reduces to the non-resonant, homogeneous granular chain, while at a very 
large mass ratio and zero initial conditions for the secondary mass the system approaches a 
model of Newton’s cradle 15 , a granular chain with quadratic onsite potential. In (Iff 


see 


also 16,17,25 ), the so-called discrete p-Schrodinger (DpS) modulation equation governing 


slowly varying small amplitude of oscillations was derived and used to prove existence of 
(and numerically compute) time-periodic traveling wave solutions and study other periodic 
solutions such as standing and traveling breathers. 

In the present setting of locally resonant granular crystals, we explore predominantly two 
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classes of solutions, namely (a) small-amplitude periodic traveling waves and (b) (dark) dis¬ 
crete breathers. To investigate these solutions at finite mass ratio (i.e., away from the above 
studied limits), we follow a similar approach and derive generalized modulation equations 
of the DpS type. We show that these equations capture small-amplitude periodic traveling 
waves of the system quite well when the mass ratio is below a critical value. We observe 
that the system admits only trivial exact bright breathers that involve linear oscillations. 
However, we use the DpS framework to prove that when the mass ratio is sufficiently large, 
the system has long-lived nontrivial bright breathers for suitable initial conditions. We also 
use solutions of the DpS equations to form initial conditions for numerical computation of 
dark breather solutions, whose stability and bifurcation structure are examined for different 
mass ratios. When the breather frequency is above the linear frequency of the resonator 
but sufficiently close to it, we identify two families of dark breathers, as is often the case in 
nonlinear lattice dynamical systems 10 . The dark breather solutions of site-centered type 
are long-lived and exhibit marginal oscillatory instability. Meanwhile, the bond-centered 
solutions exhibit real instability. This can lead to the emergence of steadily traveling dark 
breathers in the numerical simulations. In addition, we identify period-doubling bifurcations 
for the bond-centered solutions. The instability of breather solutions is also affected by the 
mass ratio. In particular, the real instability of the bond-centered breather solutions at a 
given frequency gradually becomes stronger as the mass ratio increases. 

The paper is organized as follows. Sec. [2] introduces the model, and the generalized 
DpS equations are derived in Sec. 3T In Sec. m we show that for sufficiently large mass 
ratio the equations reduce to the DpS equation derived in 15 and rigorously justify the 


validity of this equation on the long-time scale for suitable small-amplitude initial data. 
We use the modulation equations in Sec. [4] to numerically investigate small-amplitude time- 
periodic traveling waves, including their stability, the accuracy of their DpS approximation 
and the effect of the mass ratio. In Sec. [5] we show that the system admits only trivial 
exact bright breather solutions that do not involve Hertzian interactions. We then prove and 
numerically demonstrate the existence of long-lived nontrivial bright breathers at sufficiently 
large mass ratio. In Sec. [6] we construct the approximate dark breather solutions using the 
DpS equations. We use these solutions and a continuation procedure based on Newton- 
type method to compute numerically exact dark breathers and examine their stability and 
bifurcation structure in Sec. [7j Concluding remarks can be found in Sec. [8j 


2 The model 

Consider a chain of identical particles of mass mi and suppose a secondary particle of mass 
m 2 is attached to each primary one via a linear spring of stiffness K > 0 and constrained to 
move in the horizontal direction. As mentioned in the Introduction, the harmonic oscillator is 
meant to represent the primary vibration mode of a ring resonator attached to each primary 
mass or a cylindrical rod. Let u n (t) and v n (t) denote the displacements of the 77 th primary 
and secondary masses, respectively. The dynamics of the resulting locally resonant granular 
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chain is governed by 


m i 


d 2 u n 
dt 2 


m 2 


d 2 v n 
dt 2 


A(u n — i A(u n K (u n u n ), 

A'(hn - V n ). 


( 1 ) 


Here A{u n — ti n+ i)“ is the Hertzian contact interaction force between nth and (n + l)th 
particles, where (x) + = x when x > 0 and equals zero otherwise, so the particles interact 
only when they are in contact, A > 0 is the Hertzian constant, which depends on the 
material properties of the contacting particles and radius of the contact curvature, and 
a is the nonlinear exponent of the contact interaction that depends on the shape of the 
particles and the mode of contact (e.g. a = 3/2 for spherical beads and orthogonally stacked 
cylinders). Typically, we find a > 1, although settings with a < 1 have also been proposed; 
see e.g. 1261 and references therein. In writing (jl| we assume that the deformation of the 
particles in contact is confined to a sufficiently small region near the contact point and varies 
slowly enough on the time scale of interest, so that the static Hertzian law still holds 13]; this 
is known to be a well justified approximation in a variety of different settings (l, 2 . We also 
assume that dissipation and plastic deformation are negligible, which is generally a reasonable 
approximation, although dissipation effects have been argued to potentially lead to intriguing 
features in their own right, including secondary waves |27 (see also |28j|). Choosing R to be 
a characteristic length scale, for example, the radius of spherical or cylindrical particles, we 
can introduce dimensionless variables 


u n v n ~ I R a 1 A 

Un = ~R' Vn = ^R } 1 = t \j mi 

and two dimensionless parameters 

m 2 K 

P > ^ A pa-1 ’ 

mi AR a 1 

where p is the ratio of two masses and n measures the relative strength of the linear elastic 
spring. In the dimensionless variables the equations (|TJ) become 

hn (Mn-1 Un)_ j_ (Un ^n+l)-f V n ) 

pV n W) ■} 

where u n and v n are second time derivatives. In what follows, it will be sometimes convenient 
to consider ([2]) rewritten in terms of relative displacement (strain) variables x n = u n — u n -i 
and y n = v n - v n -p. 


%n 2( X n )_j_ ( X n -\-i)+ ( X n _i)_|_ K(x n Vn) , % 

, i I'j) 

PVti Vn)• 

Note that in the limit p —>• 0, the model reduces to the one for a regular (non-resonant) 
homogeneous granular chain. Meanwhile, at p —> oo and zero initial conditions for v n (t) 
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the system approaches a model of Newton’s cradle, a granular chain with quadratic onsite 

(4) 


potential, which is governed by 15,29] 


U n ~t~ (u n — i ( U n Un- 1-1 )- 


In 15 , the discrete p-Schrodinger (DpS) equation 
dA r 


2 iTn 


dr 


— (A n +1 — A n )\A n+ i — A n \ a — (A n — A n _i)\A n — A n _i 


la—1 


(5) 


has been derived at k — 1 to capture the modulation of small-amplitude nearly harmonic 
oscillations in the form 


= e{A n (r)e lt + A n (r)e lt ), r = e a H, (6) 

where e > 0 is a small parameter, A n (r) is a slowly varying amplitude of the oscillations 
and r 0 is a constant depending on a. In the next section, we follow a similar approach and 
use multiscale expansion to derive generalized modulation equations of the DpS type for (|2]) 
with finite p. 


3 Modulation equations for small amplitude waves 

3.1 Derivation of generalized DpS equations at finite p 

Using the two-timing asymptotic expansion as in [l5), we seek solutions of (]2]) in the form 

u(t) — £U(t,r), v(t) = eV(t,r) 

where r = e a ~ 1 t is the slow time, and u, v, U and V are vectors with components u n , v n , 
U m V n , respectively. The governing equations (j2]) then yield 

{d t + £ a ~ 1 d T ) 2 U = £ a ~ 1 G(U) + k(V - U) 

(dt + ^drfV = -(U-V), ^ 

P 

where the nonlinear term is given by 

G{u) n = - u n y + - (u n - u n+1 )%. 

The solution has the form 

U = u° + £ a ~ l u l + o(e Q_1 ), V = V° + e^V 1 + 

The Oth order terms satisfy a linear system, which after the elimination of secular terms 
yields 

U° = B(r) + K[A(r)e iujt + A( r )e~ iujt ], U° = B(r) --[A(r)e iut + A(r)e- iw \ (8) 

P 
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where u = \Jk + k/ p is the frequency of harmonic oscillations. This internal frequency of 
each resonator is associated with the out-of-phase motion of the displacements U and V. On 
the order 0 (e a_1 ), the system g results in 

(< d 2 + k)U 1 - kV 1 = -2 Kd T Aiue iut + c.c. + G(U°) 

d 2 V x -~{U X - V l ) = —d r Aiue** + c.c., (9) 

P P 

where c.c. denotes the complex conjugate. Let 


J(f) = 


LJ 


r» 27 r/LJ 


2W„ 


denote the projection of f(t) on e luJt and define the averaging operator as 


( 10 ) 


E(f) = 



( 11 ) 


The projection operator on all remaining Fourier modes (of the form e lJi0t , j ^ ±1, j ^ 0) is 
given by 

U h — I — E — e iwt J - e~ iujt J. (12) 

Let U\ = n h U x , V/l = U h V x . Then g yields 

(d 2 t + k)U x - nV x = n h G(U°) 
d 2 X - y h - V x ) = 0. 

Note that this equation has a unique 27r/o;-periodic solution (U^V^Y because for each j 
such that j 7 ^ ± 1 , j 0 , the matrix 


k — j 2 u 2 —k 
- n/p k/p - j 2 uj 2 


for the left hand side of the above equation associated with jth harmonic is invertible. Let 


U 1 = U\ + C°(r) + C x {r)e iujt + c.c., y 1 = V£ + D°(t) + D 1 {r)e i “ t + c.c. 
and project g on e lut , recalling that c u = y/n + k/ p: 

- ( n/p)C x - kD 1 = -2 Kd T Aioj + J(G(U 0 )) 

— (k/ p)C x — kD 1 = (2 n/p)d r Aiuj. 

This yields the compatibility condition 


2 iKtu 3 d T A = J(G(U 0 )) 


( 13 ) 


and 

D 1 = --(C 1 + 2iud r A). 
P 
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Taking the average of (j9|), we obtain 


k(C° - D°) = E(G(U 0 )) 

- ~(C° - D°) = 0. 

P 

Since p is finite, this yields C° = D°, and hence the following condition can be obtained for 
the leading order solution: 

E(G(U 0 )) = 0. (14) 


To obtain the generalized DpS equations, we now consider the conditions (13) and (14) 
in more detail. Observe that for b e M, z = re™ G C, we have 


r-2n 


E[{—b + Kze lU}t + Kze *“)“] = — / (-5 + 2 nr cost)* dt = g a (b, r), r — \z\. (15) 

2vr J o 

Here we rescaled time in the averaging integral and used the fact that the result is indepen¬ 
dent of 9 since we can always shift time when averaging. Similarly, 


r-2 7r 


J[(-b + Kze lut + Kze-**)*] = 


2tt r 


e lt (—b + 2nrcost)^_dt = zh a (b, r), r — \z\. (16) 


Defining the forward and backward shift operators 

(S + A) n = A n+ i — A n , (5 A) n = A n — A n _ i, 

we observe that 

G(U°) = -5 + (-5~U 0 )l = —5 + (—5 ~B(t) - k 5~A{r)e iu]t - A{r)e~ iuit ) a +) 


where we used the first of (J8J) to obtain the second equality. Substituting this in (13) and 
(14) and using (15), (16) with z = — 5~A and b = 5~B , we obtain the generalized DpS 
equations 

2iKuj 3 d T A = 5 + [h a (5-B,\5-A\)5-A] (17) 

and 

8 + g a {8-B, |rH|) = 0. (18) 


3.2 DpS equation at large p 


We now investigate the case of large p. Consider first the “critical” case when p = e 1 a . The 
0 th order problem is 

d?U° = k(V° - U°), d?V° = 0, 


which yields 

U° = B(t) + K[H(r)e i ^ + H(r)e-^], V° = B{t), (19) 

where we used the fact that uj = yfn when p — » oo. Meanwhile, the 0(e“ _1 ) problem becomes 

(( d 2 t +k)U 1 - kV 1 = -2iK 3/2 d T Ae iV,it + c.c. + G(U°) 
d 2 V x = t^Ae^ + c.c. 
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Note that the right hand side of the second equation has zero time average, as it should to 
be consistent with the left hand side, and the equation yields 

V 1 = -kA^Y^ + c.c. + D°(t). 

Putting this back into the first equation and projecting on , we get 


2 n 3 / 2 id r A+ n 2 A = J[G(U 0 )}, 


which is almost like the DpS equation in 15 if we set k — 1. Note, however, the additional 
term n 2 A in the left hand side and the fact that U° also includes B{r). Observe further that 
there are no conditions to determine B at this order. If B — 0, we get a (modified) DpS 
equation for A only. 

Now suppose p = £ 1-7 , 7 > a. Then the 0th order equation is the same, so the solution 
is still given by (19), while on 0(e° 1 ^ 1 ) we get 


(( d 2 +n)U 1 - kV 1 = - c l^ l2 d T Aie i ^~ Kt + c.c. + G(U°) 

d‘tV l = 0 , 


so the second equation yields V 1 = D 0 (r), while the projection of the first on e l ^ t yields 
the DpS equation for the Newton’s cradle model: 

2 K 3/2 id T A = J[G(U 0 )}- 


Note that B(t ) in (19) is again not determined at this order. Observe, however, that in the 
limit p —* oo the initial conditions r>(0) = b(0) = 0 yield v{t) = 0 (and thus V° = V 1 = 0), 
and we recover ([b]) and the DpS equation ([5]) at k = 1 : 

id T A = ou 0 5 + [\5-A\ a - 1 5-A\ (20) 


with (see |T5|) 


->a—2 


r(a?/2 + i) 


0 v^P((a + l)/2 + l)' ( } 

In Theorem 3.1 below, we justify the DpS equation ( | 20 j ) on long time scales, for suitable 
small-amplitude initial conditions. We obtain error estimates between solutions of and 
modulated profiles described by DpS. We seek solutions of (J2| and the DpS equation in the 
usual sequence spaces £ p with 1 < p < +oo. For simplicity we state Theorem 3.1 in the case 
n = 1. 


Theorem 3.1 Fix constants C r ,Ci,T > 0 and a solution A e C 2 ([0, T], £ p ) of the DpS 
equation (20). There exist constants Et > 0 and Ct > C\ such that the following holds: 

For all £ < Et and for p~ l < C^e 2 ^ 1 ^, for all initial condition (w(0), n(0), u(0), n(0)) G i p 
satisfying 

||u(0) -2eRei4(0)|| p + ||m( 0) +2elmi4(0)|| p < C j£ a , (22) 

||u(0)||p < C;E a , ||h(0)|| p < Ci£ 2 “ _1 , (23) 

the corresponding solution of satisfies for all t G [0, 

|| u(t) — 2eRe (A^e 0 ^ 1 ^ e lt )\\ p + \\u(t) + 2elm {A{E OL ^ 1 t) e lt )\\ p < C T £ a : 


\v(t)\\ p < C T £ a , \\v{t)\\ p < C T e 


2a -1 


(24) 

(25) 









Similarly to what was established in 17 for Newton’s cradle problem, Theorem 3.1 


shows that small 0(e) solutions of ([2]) are described by the DpS equation over long (but 
finite) times of order £ 1- “. However, there are important differences compared to the results 
of flT|. Firstly, the DpS approximation is not valid for all small-amplitude initial conditions, 


since one has to assume that u(0) and 7(0) are small enough (see (23)). Secondly, p must 
be large when e is small. More precisely, p must be greater than e^ l ~ a \ which scales as 
the square of the characteristic time scale of DpS. This is due to the translational invariance 
of ({ 2 ]), which introduces a Jordan block in the linearization of (2j) around the trivial state, 
inducing a quadratic growth of secular terms (see the estimate (44) below). 


Let us now prove Theorem 3.1 The main steps are Gronwall estimates to obtain solutions 
of (J2]) close to solutions of the Newton’s cradle problem Q when p is large, and the use of 
the results of 17| to approximate solutions of Q with the DpS equation. 

Equation ([ 2 ]) at k = 1 reads 


u + u — v = G(u), 

v = a (u — v ) 

where cr = 1/p is a small parameter satisfying 

a < G r £ 2(a - 1} , 


as assumed in Theorem 3.1 In addition, we have 


\\G(u)\\ p = 0(\\ur p ), \\DG(u)\\c{t p ) — 0(\\u Up 


a— 1 \ 


(26) 

(27) 


(28) 


(29) 


To simplify subsequent estimates, it is convenient to uncouple the linear parts of (]26|) and 

ariab 

R , v = R — a Q. 


(27), which can be achieved by making the change of variables 

u = Q 


The system (26)-(27) is equivalent to 


Q + Q — x G(Q + R) — u Q, 
R = x a G(Q + R), 


(30) 

(31) 


where y = (1 + a) 1 is close to unity. 


To approximate the dynamics of (30)-(31) when cr is small, we first consider the case 


cr = 0 and R = 0 of (30) (leading to the Newton’s cradle problem) and use the results of 117 


relating Newton’s cradle problem to the DpS equation. More precisely, given the solution A 
of the DpS equation considered in Theorem |3.1[ we introduce the solution Q a of 


Qa + Qa — G(Q a ) 

with initial condition Q a (0) = 2eReH(0), Q a (0) = — 2£lmH(0). According to Theorem 2.10 
of [17|, for e small enough, the solution Q a is defined on a maximal interval of existence 
(f“,f + ) containing [0, Te l ~ a ] and satisfies for all t G [0, T e 1- "] 


\Qa(t) - 2eRe(A(e a t) e l )\\ p +\\Q a (t)+2elm(A(e a t)e 


| P <Ce a . 


(32) 
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This implies in particular that 


IIQa||L°°((0,T£ 1 -“),4) + ||Qa||L° o ((0,T£ 1 - a ),^) < Ms. 


(33) 


Next, our aim is to show that Q remains close to Q a (and its DpS approximation) and R 
remains small over long times, for suitable initial conditions, e small enough and p large 


enough (i.e. a small enough). Setting Q = Q a + W in (30)-(31) yields 


W + W = N(W,R ) 

R = xvGiQ^ + W + R), 


where 


N(W, R) = -<r(Q a + x G(Q a )) + X [ G(Q a + W + R) - G(Q a ) ]-cr W. 
Moreover, from the identities 

W - Q - - (1 + o)~ l (u - v) - Qz, R = (1 + <y)~ l (era + v) 


(34) 

(35) 

(36) 


and the assumptions (22)-(23) and (28), it follows that 


11 ^(0)11,+ l|W(0)|| r <Ce“, 
l|fl(0)|| p < Cs a , J|«(0)|| p < C's 2 "- 1 . 


(37) 

(38) 


Let X = (W, R,W, R) and ||X|| P denote the sum of the £ p norms of each component. We 
shall now use Gronwall estimates to bound ||X(t)|| p on the time scale considered in Theorem 

EU 


The solutions of (34)-(35) corresponding to initial conditions satisfying (3T)-(38) are 


defined on a maximal interval of existence (t m in,^max) with t min < 0 and f max < t + , which 


depends a priori on the initial condition and parameters (in particular, e). From (37)-(38), 
one can infer that ||X(0)|| P < Me (the size of Q a in (33)) for e small enough. Let 


t E = sup {t G (0, min(T e Q ,t max )): Vs e (0,t), ||X(s)|| p < Me}. 


From (33) and the triangle inequality, we obtain 

\\Q*(t) + W(t) + R(t)\\ p <2Me, Vf e [0,fj. 


In addition, from definition (39) we have either 


L: < min(Ts 1 ",i m „) and ||A'(i e )|L = Me 


(39) 


(40) 


(41) 


or 


t F = T £ 


1 —a 


< t 


max 
1 — a\ 


(42) 


(if ||X|| p is bounded on [0, f max ), then t max = t + > Te 1 “). Integrating (35) twice yields 


R(t) = r(0) t + R( 0) + x a / (t — s) G(Q a + W + R)(s) ds. 

Jo 


(43) 
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Using the fact that t e <Te 1 ", estimate (40) and the first bound in (29), one obtains from 
the above identity the following inequality: 

VtG [0,4], ||-R(£)||p < T £ 1_Q ||i?(0)||p + ||.R(0)||p + Ciae 2 ~ a . (44) 

(45) 

(46) 

(47) 


Then the assumption (28) and the property (|38|) yield 
Similarly, we have 


Vte [0,y, \\R(t)\\ p <Ce a . 


R(t ) = R{ 0) + ya / G(Q a + W + R)(s) ds, 

Jo 

which implies 

Vte[0,t e ], \\R(t)\\ p <\\R(0)\\p + O(ae)<C£ 2a - 1 . 
Moreover, using Duhamel’s formula in (p4|) yields 


W{t) = cost W(0) + sin tfU(O) + / sin (t - s) [N(W, R)](s) ds. 


Recalling that t e < Te 1 ", using (37), and using (33), (28), (29) and (45) to estimate 
N(W,R) from the definition (36), we get 


Vte[0,t e ], \\W(t)\\ p <C 1 £ a + C 2 e a - 1 / ||VU(s)|| p ds. 

Jo 

By Gronwall’s lemma we then have 

Vte [0,f £ ], \\W(t)\\ p <C 1 e ChT e a . 

Similarly, we obtain 

W(t) = -sint W(0) + cost 1U(0) + [ cos (t - s) [N(W, R)}(s) ds. 


(48) 


(49) 


Using the same estimates as the ones involved in proving (48) and estimate (49), one can 
show that the above identity yields 


Vf6[0,iJ, \\W(t)\\ p <Ce° 


(50) 


Summing estimates (45), (47), (49), ( J50j) , we find that ||X(i)|| p = 0(e Q ) < Me for all 
t G [0, t £ ) if e is small enough. Consequently, the property pi) is not satisfied, which implies 
that (|42l) must hold instead. We therefore have 


t E — T e 


1—a 


in estimates (45), (47), (49) and (50). Combining these estimates with the DpS error bound 
(32) and the assumption (28), we deduce the error bounds (24), (25) from the identities 

u(t) - 2e Re (A(e a ~H) e u ) = W(t) + R(t) + Q a (t) - 2e Re (A(e a ~H) e u ), 

u(t ) + 2e Im (^(e" -1 ^) e lt ) — W(t) + R(t ) + Q a (t) + 2e Im (^(e" -1 ^) e lt ), 

v = R — a Q a — cr W. 


This completes the proof of Theorem 3.1 
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4 Time-periodic traveling waves 

We now use the DpS equations to investigate time-periodic traveling wave solutions of our 
system. In what follows, it will be convenient to use the strain formulation ([3]) and also 
rewrite the DpS equation (17) in terms of strain-like variables 


2 iKu> 3 d T 8 A n = 5 + [h a (5 B n ,\S A n \)5 A n ] - 5 + [h a (5 B n _ u \8 A n _i|)<5 Ai-i]- (51) 

A special class of solution of (51), (18) has the form of a periodic traveling wave 

5~A n (r) = S~B n (r ) = b , (52) 

where the frequency fl depends on the amplitude a > 0, the constant b and the wave number 
q through the nonlinear dispersion relation 


Q = 




• <? 
sm 2 


h a (b, a). 


(53) 


Note that definition of h a in (16) requires that b < 2 na. In particular, the dispersion relation 


has a closed form when 6 = 0. In fact, as shown in (151 for k— 1, in this case 
h a (0,r) = 2 a ~ 1 c Q K a r a - 1 , 
and the dispersion relation is given by 


ir(l/2)r(a/2 + l) 
Ca- 7 T r((a + l)/2 + 1) ’ 


fin — 


c a 2 a (na) 


a—1 


( 


. q 

sm - 
2 


Equations (|8l), (52) and (53) yield the following first-order approximation of the periodic 


traveling wave solutions of system (j3]): 

xt n (I) = zb + 2nea cos (nq — u> tw t + 0), y(t) = eb 


2 nea 
P 


cos (nq — u tw t + (j>), (54) 


where oo tw is tl ie traveling wave frequency given by 


CUt w — fin' 0 -)- LO — to -)- 


2e' 


a—1 




• q 

sm - 
2 


h a (b, a) 


(55) 


To investigate how well the dynamics governed by the DpS equations approximates the 
traveling wave solutions of (J3|, we consider initial conditions determined from the first-order 
approximation 


x n PP (t) = B n {e a 1 t) + keS A n (e a l t)e luJt + c.c. 

y* pp {t) = E5~B n {e a -H ) - —5~A n {E a -H)e iuit + c.c. 

P 

at t — 0, along with initial velocities given by 

x n (0) = £{e ci ^ 1 S~ B n ( 0) + KE a ~ l 5~ A n ( 0) + inu;8~ A n (0) + c.c} 


(56) 


K 


. K 


y n (0) = £{E a 5 B n ( 0)- £ a A n (0) - i-u8 A n (0) + c.c}. 

P P 


(57) 
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Here S A n (0) are set to be small perturbations of (52) at r = 0 and 0 = 0: 

(TH n (0) = a[(l + Ci 1} ) cos (nq) - i{ 1 + C i 2) ) sin (nq)], 


( 58 ) 


where Cn'^ and (i 2 ' 1 are uniformly distributed random variables in [—(, £] with small ( > 0. 
Let C = g a (b,\a\) > 0 for given constants a and b < 2/c|a|. Then 5~~B n (r) is determined 


from (18) by numerically solving 


g a {S-B n {r),\S-A n {r)\) = C. 


(59) 


From the definition of g a in (15), it is clear that the function b t—)■ g a (b,r ) is decreasing on 


(—oo,2ftr) and satisfies lim^.oo g a (b, r) = +oo and g a (b,r) = 0 for b > 2 nr. Consequently 
equation (59) admits a unique solution 5~B n (r) G (—oo, 2n\5~A n (T)\). In particular, we 


obtain the value of 5 B n { 0) by solving (59) at r = 0. Note that S A n (0) and S B n { 0) can 


be computed exactly using DpS equations (51), (59), although their contribution is negligible 


since they correspond to higher order terms in the expression of initial velocities x n (0) and 
y n ( 0). In particular, when the perturbation is zero (i.e. ( = 0), the initial condition simplifies 
to 


x n (0) = eb + 2ehiacos(nq + 0), y n ( 0) = eb — 


2 EKCL 


cos (nq + 0), 


P 

i n (0) = 2u) tw K,ea sin(ng + 0), yJ 0) = — s j n ( n q + 0). 

P 


(60) 


Integrating (|3]) numerically on a finite chain with these initial conditions and using periodic 
boundary conditions x 0 = x^, Xjv+i = Xi, we can compare the solution of (J3]) (referred to in 
what follows as the numerical solution) with the ansatz (54) when ( = 0 or the ansatz (56) 


when C > 0. The latter is obtained by solving (51), (59) with periodic boundary conditions 


5~A 0 (r) = 5 ~An(t), <5 - Htv+i(t) = 5~Ai(t) 1 using the Runge-Kutta method and a standard 
numerical root finding routine to determine 5~B n {r) from (59) for given h _ H n (r) at each 
step. 


4.1 Numerical traveling waves and the DpS approximation 

We now consider a locally resonant chain with N — 50 masses. In what follows, we set 
k — 1, noting that other values of this parameter can be recovered by the appropriate 
rescaling of time and amplitude. To investigate the accuracy of the DpS approximation 
of the small-amplitude traveling waves, we first fix mass ratio p = 1/3 and normalize the 
traveling wave solution (52) by fixing a = 1, 0 = 0 and 6 = 1. The linear frequency is given 
by c o = k + n/p = 2. 

In the first numerical run, we set q — 7 t/5 , £ = 0 and consider the traveling wave 
with frequency u> tw = oo + 0.001, so that £ ~ 0.057 from the nonlinear dispersion relation 
(55). Equation (54) then yields the amplitude of x}'“ approximately equal to 0.114, which 


corresponds to the small-amplitude regime. As shown in the left panel of Fig. [TJ the agree¬ 
ment between the numerical and approximate solutions is excellent, even after a long time 
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t = 100 T to , where T tw = 2n/u tw = 3.126 is the period of the traveling wave. The relative 
errors of the approximate solutions 

E x(t) = 2 ~H a £ SP W ^ x n(t)\\oo and E y {t) = ^—\\y“ pp (t) - y^t)^ (61) 

are less than 8% at the final time of computation and remain bounded throughout the 
reported time evolution, as shown in the right panel of Fig. [lj In the second numerical run, 




Figure 1: Left plot: strain profiles of small-amplitude approximate solution (54) (connected stars) and 
numerical solution of © (connected squares) at t = 100 T tw « 314. Right plot: the relative errors E x {t) 
(black curve) and E y (t) (grey curve) of the DpS approximation. Here (j> = 0, k = 1, q = 7r/5, b = 1, a = 1, 
C = 0 and uJt w = w + 0.001. 


we set C — 0.01 for the perturbation in ( |58| ), while the other parameters remain the same. 
The agreement between the numerical and approximate solutions is still excellent over the 
same time interval (see the left plot of Fig. [2]). Moreover, the right plot of Fig. [2] shows the 
normalized differences 

Exit) = ^ a W x n d ( t ) -Zn(*)llo° and Ev(t) = 2 ^\\y£ d (t)-y n (t)\\oo (62) 

between a perturbed traveling wave (x^ d , y^ d ) (numerical solution obtained for the perturbed 
initial condition) and the unperturbed numerical traveling wave solution (x n , y n ) shown in 
Fig.0 with the wave number q = 7t/5. As shown by Fig. [2j the initial perturbation is 
not amplified at the early stage of the numerical integration of ([3]) (for t < 50 ~ 16 T tw ). 
However, the subsequent growth of perturbations indicates the instability of the traveling 
wave solution. 

In the next computation, we increase the wave number up to q = 47 t/5 and keep all 
the other parameters the same as before. Now the asymptotic scale becomes quite small as 
e ~ 6.37 x 10“ 4 , which yields a very small amplitude of x l ™ ps 0.0013. As shown in Fig. [3j 
the DpS equations can successfully capture the dynamics of the small-amplitude traveling 
wave solution of the original system. In addition, the small-amplitude traveling wave with 
q = 47 t/5 appears to be stable on the interval t G [0, 600]. However, this result may be linked 
with the very small traveling wave amplitude, and instabilities might appear on longer time 
scales. 

It is worth pointing out that the time scale of the validity of the DpS approximation 
depends not only on the asymptotic scale e but also on the wave number q. To illustrate 
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Figure 2: Left plot: strain profiles of small-amplitude approximate solution (connected stars) from the 
ansatz ( |56| and perturbed traveling wave solution (connected squares) at t = 100 T tw ss 314. Right plot: the 
normalized differences E x (t) (black curve) and E y (t) (grey curve) of the perturbed traveling wave and the 
unperturbed solution shown in Fig. [l] Here £ = 0.01 and all the other parameters are the same as in Fig. [l] 
A growth of the perturbations can be clearly observed in the dynamics. 



Figure 3: Top panels: strain profiles of small-amplitude approximate solution (connected stars) from the 
ansatz ( |56[ ) and numerical results (connected squares) with both unperturbed (left) and perturbed initial 
conditions with ( = 0.1 (right) at t = 200 T tw « 628, respectively. Bottom panels: left plot shows the 
relative errors E x (t) (black curve) and E y (t) (grey curve) of the DpS approximation. Right plot shows the 
normalized differences E x (t) (black curve) and E y (t ) (grey curve) between the perturbed and unperturbed 
traveling waves. Here q = 47r/5 and all the other parameters are the same as in Fig. |T] 
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this, we first consider the wave number q = 7t/ 5 but increase the traveling wave frequency 
up to uJt w = uj + 0.003, which yields e ~ 0.514 and the amplitude of around 1.028. As 
revealed by the left plot in Fig. |dj the DpS equations fail to describe the dynamics of (|3]) 
appropriately soon after we start the integration. We then increase the wave number up 
to q — 47t/ 5 but choose cu tw = c o + 0.009 yielding e ~ 0.0516, which is even smaller than 
the asymptotic scale e used in Fig. [lj However, notable difference between the traveling 
wave patterns of the numerical and approximate solutions are observed over the same time 
interval [0,314] (see the right plot of Fig. [I]). 



Figure 4: Left panel: results of the simulations with the same parameters as in Fig. [l] except for ui tw = 
uj + 0.003 and the snapshot is taken at t = 10 T tw k, 31.4. Right panel: results of the simulations with the 
same parameters as in the left panels of Fig. [3] except for u! tw = w + 0.009 and the snapshot is taken at 
t = 314. 


4.2 Effect of mass ratio p 


To investigate the effect of the mass ratio p on the validity of the DpS approximation of the 
small-amplitude traveling waves, we fix q = 7t/5, e = 0.01 and keep the other parameters 
the same as before. We choose mass ratio p = 3 and the traveling wave frequency in (55) 
is now given by u tw cs uj + 0.002, where the linear frequency is uj — 1.1547. The results 
of the simulations are shown in Fig. [5} We observe again that the agreement between the 
numerical and approximate solutions remains excellent over the time interval [0,100 T tw \, and 
the numerical traveling waves appear to be stable over the time interval [0, 100]. However, 
we note the growing trend of the difference between the perturbed and unperturbed traveling 
wave solutions after t ta 100, which illustrates the instability of the traveling wave. 

We further increase p to the critical value p c = e l ~ a = 10. The linear and traveling 
wave frequencies are given by uj = 1.0488 and ujt w = 1.0517, respectively. As shown by the 
left panel of Fig. [ 6 j the wave form of y n is no longer sinusoidal at t — 100 T tw , while the 
wave form of x n remains sinusoidal and matches its DpS approximation very well over the 
time interval [0,100 T tu: ]. We observe also a growing trend in the relative error between the 
numerical and approximate solutions, despite the structural similarities of the profiles. To 
investigate the case when p is large ( 0 (£ 1-7 ), 7 > a), we now set p = 1000 and keep all 
other parameters the same as before. The linear frequency is now given by uj — 1.0005 and 
approximate frequency is oj tw = 1.0038. Again, over the integration time interval [0,100 T tw ] 
the wave form of x n remains sinusoidal and close to the ansatz (54) but the waveform of 
y n is non-sinusoidal and significantly deviates from the DpS traveling wave approximation 
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Figure 5: Results of the simulations with the same parameters as in Fig. [ 3 ] except for p = 3, q = 7 t/5 
and e = 0.01. The strain profiles of both panels correspond to snapshots taken at t = 100 T tu , and the 
perturbation added in the initial condition is £ = 0.01. 




Figure 6: Results of the simulations with the same parameters as in the left panel of Fig. [H] except for 
p = 10 (left panel) and p = 1000 (right panel). The strain profiles of both panels correspond to the snapshot 
taken at t = 100 T tw . The insets in the bottom plots show the enlarged plots of the relative error E x (t). 
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(see the right plots in Fig. [6]). These results are consistent with the discussion in Sec. 3.1 


since the derivation of the DpS equations (17)-(18) requires that p is below the critical value 
p c = e 1 '" for given e. Note that if we further set 6 = 0, the DpS equation ( flr[ ) reduces 
to the one capturing the traveling waves in Newton’s cradle problem in [15]. Numerical 
results (not reported here) reveal that the numerical solution of x n is sinusoidal and very 
well approximated by the traveling wave a n sat z (54). However, the difference between exact 
and approximate solutions of y n is very substantial. Here the structural characteristics of 
the solution for y n are no longer properly captured by the DpS approximation. 

These numerical simulations reveal that the validity of the DpS approximation at fixed 
wave number q is very sensitive to the mass ratio p. When p is relatively small, the generalized 
DpS equations ([I7|)-(fl8| successfully capture the dynamics of periodic traveling waves. When 


P > Pc = £ 


1—a 


an increasing deviation between the exact and approximate solutions of y r 


emerges. However, in all cases, the agreement of approximate and numerical solutions for 
x n remains excellent over a finite time interval. In addition, traveling wave instabilities can 
be observed depending on the values of p, wave number q, wave amplitude and time scales 
considered. 


5 Bright Breathers 


Bright breathers are time-periodic solutions of (|2]) which converge to constants (zero strain) 
at infinity, i.e. 

u n ,v n —* c±, as n —> ±oo (63) 


uniformly in time. In this section we examine the existence of either exact or long-lived 
bright breather solutions of (|2j) . The second class of solutions refers to spatially localized 
solutions of ([2]) remaining close to a time-periodic oscillation over long times. 

We begin by noting the existence of trivial exact bright breather solutions of ([2]) for which 
particles do not interact, i.e. (u n _i(f) — w n (t))“ = 0 for all t and n. This is equivalent to 
having 

u n ~i(t) < u n {t ) Vt G M, Vn G Z, (64) 


U n — t^iVn 
pV n ^-(®n 


The time-periodic solutions of (65) read 


u n ), 

V n ). 


(65) 


Unit) = 


P 


1 + P 


COS {ut + 4>n) + b n , v n {t) = 


1 + P 


a n cos {ut + 4> n ) + b n , 


( 66 ) 


where we can fix a n > 0 and denote oj 2 = k (1 + 1/p). Bright breather profiles are obtained 
for 


lim a n = 0, lim b n — c±. 

n —^ioo n—>d=oo 


(67) 


A solution of (J2]) is obtained if and only if the constraint (64) is satisfied, which is equivalent 
to 


6n-i > 


P 


1 + P 


[a 2 n + - 2a n a n -i cos (0 n - 0 n _i)] 1/2 Vn e Z. 


( 68 ) 
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For ( a n ) n£ z G ^i(Z), this is equivalent to assuming 


bn 


d n + 


P 


1 + P 


n 

E 

k=—oo 


Wl + al-i- 2a fc a fc _i cos (0 fc - 0 fc -i) ] 1/2 , 


where (d n ) n£ z is a nondecreasing sequence converging as n —> ±oo. It is clear from this 
expression that (b n ) n& i corresponds to a kink prohle. Moreover, fixing — (j) n -\ = 7r, we 
can simplify the above expression to obtain 


b 


n 


d n + 


P 

1 + P 


n —1 

(o n + 2 ''y ^ a . 

k=—oo 


We now prove the following. 


Theorem 5.1 The trivial bright breather solutions defined by (66)-(61)-(68) are the only 
time-periodic bright breather solutions of 


To prove Theorem |5.1| we follow the method of the proof of nonexistence of breathers in 
FPU chains with repulsive interactions given in [l6|. Suppose (u n ,v n ) is a bright breather, 
i.e. a T-periodic solution of (j2]) satisfying (63). Adding the equations in (J2]), one can see 
that the bright breather solution must satisfy 


hn + pV-n (pin— 1 ^ri )(+n ^n+ 1)-|-' 


Integrating (69) over one period, we obtain 

+n+l +n; 


F n = — / (u n _i(t) - u n (t))°dt. 


(69) 


(70) 


Note that F = F n is independent of n and one can show that it vanishes. Indeed, by (63), 
we have 


lim || u n — u n — \ | — 0. 

n^-=t<xi 


(71) 


Meanwhile, 


F\ = — I (t) - u n (t))“dt < ||u„_i(i) - ii„(4)||S»( 0i t) 


(72) 


for all n. Taking the limit n —)■ ±oo in (72) one obtains F = 0. Consequently, for each n, 
we have 


{u n -i(t) - u n (t))°fdt = 0 


(73) 


and since F n = ( m„_i — «„)" is non-negative, continuous and T-periodic, we have F n = 0 for 
all t and n. Thus (u n ,v n ) satisfies the linear system (65) and corresponds to a trivial bright 
breather solution. This completes the proof of Theorem 5.1| 


19 
















In what follows we show that, although nontrivial time-periodic bright breathers do not 
exist for system ([2]), long-lived small-amplitude bright breather solutions can be found when 
p is large. This is due to the connection between ([ 2 ]) and the DpS equation (20) established 
in section 3.2 Equation (20) admits time-periodic solutions of the form 


A,(T)=£a„e‘"' , W"-y 


(74) 


where a = (a n ) neZ is a real sequence and £ G M an arbitrary constant, if and only if a satisfies 

(75) 


= -a + [|ra|“-‘ra]. 


In particular, nontrivial solutions of (75) satisfying lim^ioc a n = 0 correspond to bright 


breather solutions of (20) given by (74). These solutions are doubly exponentially decaying, 


so that they belong to i v for all p G [1, 00 ]. They have been studied in a number of works 
(see 25 and references therein). The following existence theorem for spatially symmetric 


breathers has been proved in 1251 using a reformulation of (75) as a two-dimensional mapping. 


Theorem 5.2 The stationary DpS equation (75) admits solutions a l n (i = 1,2) satisfying 


lim a* =0, 

ra->±oo n 


(-l) n <>0, 


Kl > K-il 

a 2 = —a 2 _ 


for all n < 0, 


and a l n = a)_ n , a z n = — a(_ n+1 , for all n G Z. 

Furthermore, for all (3 G (0,1), there exists n 0 G N such that the above-mentioned solutions 
a l n satisfy, for i — 1,2: 

Vn > no, \a l n \ < (3 1+a " "°. 


Considering the bright breather solutions of (20) given by (74) with a = a 1 , e = 1, and 
applying Theorem 3.1, one obtains stable exact solutions of equations (|2j) , close to the bright 


breathers, over the corresponding time scales. This yields the following result formulated for 
k = 1 . 

Theorem 5.3 Fix constants C r ,C[,T > 0. Consider a solution C = (<) nez (1 = 1 , 2 ) of 


the stationary DpS equation (75) described in Theorem 5.2. There exist £t, Ct > 0 such that 
for all £ G (0, £t] and for p _1 < for all initial condition of (|2]) in C p satisfying 

||«(0) - 2ea l \\ p + ||m(0)|| p + ||v(0)|| p < Ci£ a , ||u(0)|| p < Ci£ 2a-1 , (76) 

the solution to equation (|2]) is defined at least for t G [O,^ 1 ^"] and satisfies 

|| u(t) — 2 e a 1 cos (fl 6 f)||p + || u(t) + 2 sod sin (h2 fc f)|| p < CtE 0, , for all t G [0, Te 1 ^ Q ], (77) 


with D b = 1 +uj 0 £ a , w 0 defined in (21) and 


< C T e a , 


< C T E 2a -\ for all t G [0, Te 1 ^}. 


(78) 
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It is important to stress the differences between the long-lived bright breather solutions 


provided by Theorem |5.3| and the trivial bright breathers analyzed at the beginning of this 
section. The oscillations described in Theorem 15.31 are nontrivial in the sense that Hertzian 
interactions do not vanish identically. In addition, they are truly localized (for p < oo) 
whereas the trivial exact breathers are superposed on a nonvanishing kink component b n . 
Moreover, the (approximate) frequency of long-lived bright breathers satisfies 0 < flj,—1 = 
0(e a ~ l ) for breathers with amplitude 0(e). For trivial exact bright breathers, the frequency 
uj is independent of amp litud e and satisfies 0 < uj — 1 = 0(1/p) when p is large. Under the 
assumptions of Theorem 5.3, we have 1/p = 0(£ 2 (" -1 )) and thus Clb > u> for e small enough. 


We now investigate the behavior of long-lived bright breathers numerically. Fixing T = 
2h/ujq = 19.4158, C r — C t — 1 and choosing e = 0.01, we consider a locally resonant chain 
of N = 40 masses with the mass ratio p = 1000 so that the inequality p~ x < C r £ 2< ' a ~ 1 ' > in 
Theorem 5.3 is satisfied. We then integrate the system ([2]) over a long time interval [0, 80 T ], 
starting with initial conditions 


w(0) = 2 ea\ u(0) = 0, u(0) = t)(0) = 0 


where a 1 is the numerical solution of 
simulation yields spatially localized solutions of 
oscillation 


(79) 


The numerical 


obtained using the method in 16 

([2]) that stay close to the time-periodic 


u(t) = 2ea l cos(Hfct), u(t) = — 2ea l sin(Oftt), v(t) = v(t) = 0 (80) 


with i — 2 over the time interval [0, Te l ~ a ] = [0,194.158], as can be seen in the inset of the 
right panel of Fig [TJ Note that the comparison is made at times corresponding to multiples 
of Tf, = 2Ti/Qb = 6.0862. To measure the relative difference of the numerical solution of (|2]) 
and the time-periodic oscillations (|80|, we define the rescaled ^-norms as follows: 


E u (t) = 


|u(t) -U(t)||] 


Eu(t) = 


K*) -“Wlh 


(81) 


and 


E v (t) = 


|v(i) ~v(t)\\i 


Ey(t ) = 


M*) -£(*)! b 

(-2a—1 


(82) 


The fact that those rescaled norms remain small over time interval [O,!^ 1 “] is consistent 


with the result of Theorem 5^ and confirms the existence of the long-lived bright breathers. 
However, at larger time, part of the energy is radiated away from the vicinity of the initially 
excited sites. As a result, we observe the breakdown of the localized structure for a long¬ 
time (t Te l ~ a ) evolution, and the solution profile spreads out and eventually approaches 
a kink-type structure shown by circles in the left panel of Fig [TJ This is associated with the 
growing magnitude of v n during the simulation. 


6 Approximate dark breather solutions 

We now turn to dark breather solutions, which, as we will see, are fundamentally different 
from the waveforms considered in the previous section and are not excluded by the results 
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Figure 7: Left panels: snapshot of numerical solution o f (pj ) at t = 3071 < 10 T (blue stars) and t = 80 T 
(red circles), respectively, starting with initial conditions (79) with i — 2 (black squares). Right panel: time 
evolution of the rescaled £i-norms defined in (81) and where E u (t) is represented by connected dots 
(black curve), Ey,(t) by connected stars (green curve), E v (t) by connected circles (red curve) and Ey(t) by 
connected pluses (blue curve). The inset in the right plot shows time evolution of the same fi-norms when 
t < Te x ~ a . 


of Theorem 5.1 To construct approximate dark breather solutions, we start by considering 
standing wave solutions of the generalized DpS equations ([T7|)-([l8| in the form 


5-A n (T) = 5-a n e i ^ T+<t>) ( 


CLn, £ 


5 B n (r) = 5 b n , 


(83) 


where a n and b n are time-independent. Introducing LUb = to + Qe a we find that the 
first-order approximate solution of ([2]) reads 


u s ™{t) = eb n + 2 nea n cos(u)bt + </>), 

, 2k 

v™(t) = eb n - £a n cos (oj b t + 0). 

P 

Substituting (83) in the generalized DpS equations (|I7|)-(18) we obtain 

P&n ha^S b n+ i, 1 5 ®n +1 h a (5 b n , |5 a n 

5 + g a (5~b n , \8~a n \ ) = 0, 

where /i = 2 kuj 3 Q = 2K,u> 3 (uJb — wk 1 ^. 

Following the approach in 115], for p ^ 0 one can further show that a n = 
b n = \p\ T ^b n satisfy the renormalized equation 

- sign(p)a n = h a (S~b n+1 , \8~a n+1 \)S~a n+1 - h a (8~b n , \8~a n \)8~a n 
8 + g a (8~b n , \5~a n \) = 0. 


(84) 


(85) 


( 86 ) 


where sign(p) = 1 for ji > 0 and sign(/i) = — 1 for /i < 0. For simplicity we drop the tilde 


in (86) in what follows. Numerical results suggest that a nontrivial solution for {<5 a n } can 


be found if and only if p > 0. Thus it suffices to consider the case p = 1. It is convenient 


to rewrite (86) in terms of 8 a n and 8 b n by subtracting from the first equation in (86) the 
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same equation at n — 1: 


5 d n h a (S fe n _|_i ? |5 u n _|_i|)5 u n _|_x 2 h a (5 b n5 15 u n |)5 u n -\- h a (5 6 n _x, |5 &n—1|)^ &n—i 

5 r a(^ ^n+l> ^n+l|) Qol^ ^n? ^n|)’ 

(87) 


We now use Newton’s iteration to solve (87) numerically for 5 a n , 5 b n , n = —N ,..., N, 
with periodic boundary conditions (<5 _ a_N_i = 5~a m and 5~aN + i = 5~a_N)- Note that the 


associated Jacobian matrix is singular due to the structure of second equation in (87), and 


therefore an additional constraint is necessary. It is sufficient to fix S b -n = c, where c is a 
constant. Since we are looking for dark breathers, it is natural to consider initial values of 


~a n in the form 


5 a° = (—l)"tanh(n — n 0 ) 


( 88 ) 


where no is an arbitrary constant corresponding to spatial translation. One can then use 
the second equation in (87) to solve for initial guess of 5~b^, n = — N + 1,..., N. A standard 
Newton’s iteration procedure of the system (87) with 4IV + 1 variables S~b-N + i ,..., 5~b]y, 


<5 a_ at, ... ,5 ajv is then performed with the tolerance of 10 . Setting n 0 = 0 in (88) results 


in a site-centered solution shown in the left panel of Fig. [8j whereas the bond-centered 
solution corresponds to no = 1/2 shown the right panel. Typical breather waveforms of both 
bright 


12 and dark 18 type come in these two broad families 10 . 
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Figure 8: 

solution. 


Left panel: A site-centered solution of (87) with S b_N 


c = 0. Right panel: bond-centered 


Once the Newton’s iteration converges to a fixed point {(J a*, S &*)}, we can compute 
the first-order approximate solution of ([3]) given by 

x s n(t) = \fi\^ {$~b* n + 2 k 5~ a* n cos(uj b t + (p)}, 

_j_ 2 k _ (89) 

3C , ( i ) = H“ -1 U K - $ a* n cos(u b t + 0)}, 


where fi = 2 K 0 J d '{uj b — u>). 
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7 Numerically exact dark breathers and linear stability 
analysis 


Having constructed the initial seed (89), we can compute the numerically exact dark breather 


solution of system ([3]) with periodic boundary conditions. Let x(t), y(t ), x(t ) and y(t) denote 
the row vectors with component x n (t), y n (t), x n (t ) and y n (t), respectively. Let Z(t ) = 
(x(t),y(t)). We seek time-periodic solutions (Z(t),Z(t)) of the Hamiltonian system (|3]) 
satisfying the initial conditions (Z(0), Z( 0)). For a fixed period of the dark breather solution 
given by T |\ = 2tt/ ujb, where u is the breather frequency, the problem is equivalent to finding 
the fixed points of the corresponding Poincare map P Tb [(Z(0), Z(0)) T ] = (Z (T b ) , Z(T b )) T . 
Since the system (J3]) has the time-reversal symmetry, we can further restrict the solution 
space by setting Z( 0) = 0. 

We use a Newton-type algorithm (see, for example, Algorithm 2 in (30)) to compute 
the fixed point. More precisely, let AZ(0) be the small increment of the initial data that 
needs to be determined. It is then sufficient to minimize \\PT b [(Z(9) + AZ(0), 0) T ] — (Z(0) + 
AZ(0),0) / || at each iteration step. Notice that Pr 6 [(kF(0), 0) T ] can be approximated by 
PT b [(Z( 0), 0) T ] + M(Tb)(AZ(0), 0) 1 for sufficiently small AZ(0). Here M[t) is the associated 
monodromy matrix of the variational equations satisfying 


= J(Z(t),Z(t))M(t), M( 0) = /. 


(90) 


where J(Z(t ), Z{t )) is the Jacobian matrix of the nonlinear system ([3j) at ( Z(t ), Z(t)), and I 
is the identity matrix. The Jacobian for the Newton’s iteration is then given by I—M (7J) and 
it is singular since it can be shown that M{Tb ) has an eigenvalue pair equal to 1. To remove 
the singularity, we impose the additional constraint that the time average of X\{t) + Axi(t), 
the first component Z(t) + AZ(t), is fixed to be (approximately) zero. Observing that 
(Z(t) + A Z(t), Z(t ) + A Z(t)) T « (Z(t), Z(t)) T + M(t)(Z( 0), 0) T , we obtain 


wr / Xi(t)dt+^f Mi(t) ■ AZ(0)dt — 0, 
1 b Jo J o 


(91) 


where M\(t ) is the first row of M(t). 

In the results discussed below we set k — 1. To characterize the solution, we define the 
vertical centers of the solution for x and y components fl8l, 


a = 


su Pte[o,T 6 ] xi(t) + inf t 6 [ o,T 6 ] X!(t) 


Cy = 


suPteio.T,] Vi(t) + inf tG [o,T 6 ] Vi (t) 


(92) 


and the amplitudes of the breather, 

R = su Ptg[o,r 6 ] X!(t) - inf fg [ 0|Tfc ] xi(t) 


„ _ su Pte[o,T b ] Vi(t) - infte[o,T 6 ] Vi(t) , . 

I\y . W'jJ 


Note that the vertical center C x is approximately zero due to the constraint (91). However, 


one can fix any other value of C x by replacing zero in the right hand side of (91) by C x . To 
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further investigate the long-term behavior of the dark breather solution, we introduce the 
relative error 


E b (t) = II Z(mT b ) - zm^/wzmu (94) 


where m = \t/T b \ and Z(mT b ) = (x(mT b ),y{mT b )) represents the strain profile after inte¬ 
grating ([3j) over m multiple of time periods, starting with the static dark breather Z( 0) as 
the initial condition. 

We first consider the mass ratio p = 1/3, so that the linear frequency is tv = 2. We 
choose a value of tv b that is slightly greater than this value but close enough to it in order to 
obtain a good initial seed with a small amplitude. Once the Newton-type solver converges to 
a numerically exact dark breather solution, we use the method of continuation to obtain an 
entire family of dark breathers that corresponds to different values of tv b . Sample profiles of 
both bond-centered and site-centered dark breather solutions with tv b = 2.05 along with the 
DpS approximate solutions (89) are shown in Fig. 19} The amplitudes K x and K y increase 


with co b , and the solution approaches zero as tv b —> tv. 



Figure 9: Left panel: a bond-centered solution of dark breather solution (connected stars) with frequency 
u> b = 2.05. Connected squares represent strain profile after integration over T b and circles represent the DpS 


approximate solution from the ansatz (89). Right panel: site-centered solution. The relative errors E b (T b ) 
of both site-centered and bond-centered solutions are less than 4.5 x 10~ 9 . Here k = 1 and p = 1/3. 

The linear stability of each obtained dark breather solution is examined via a standard 
Floquet analysis. The eigenvalues (Floquet multipliers) of the associated monodromy ma¬ 
trix M(T b ) determine the linear stability of the breather solution. The moduli of Floquet 
multipliers for the site-centered and bond-centered solutions of various frequencies are shown 
in Fig. [TOj along with the numerically computed Floquet spectrum that corresponds to the 
sample breather profile at iv b = 2.09. If any of these Floquet multipliers A* satisfies A/ > 1, 
the corresponding breather is linearly unstable. We observed two types of instabilities in 
this Hamiltonian system. The first one is the real instability , which corresponds to a real 
Floquet multiplier with magnitude greater than one; an example is shown in the right top 


plot in Fig. 10 for the bond-centered dark breather with tv b = 2.09. The second type is the 
oscillatory instability, which corresponds to a quartet of Floquet multipliers outside the unit 


circle with non-zero imaginary parts (see the right bottom plot of Fig. 10 for the site-centered 
dark breather with tv b = 2.09). 

Numerical results reveal that the site-centered dark breathers appear to exhibit only 
oscillatory instabilities. These marginally unstable modes emerge at the beginning of the 
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Figure 10: Left plots: moduli of the Floquet multipliers versus frequency Ub for the bond-centered (top) 
and site-centered (bottom) breathers. The Floquet multipliers for co b = 2.09 in the complex plane are shown 
in the respective plots on the right. Here n = 1 and p = 1/3. 


continuation procedure but remain weak until u b reaches uj b ~ 2.078. As shown in the right 
plot of Fig. 11 the relative error 77 b (1200) stays below 7 x 10~' when oj b < 2.077 but increases 
dramatically afterwards. The results are consistent with the Floquet spectrum shown in the 
left plot of Fig. 11 In fact, beyond the critical point of uj b ~ 2.078, we observe the emergence 
of a new and stronger oscillatory instability mode that corresponds to two pairs of Floquet 
multipliers distributed symmetrically outside the unit circle around —1 (see also the right 
bottom plot of Fig. 10 for uj b = 2.09). Representative space-time evolution diagrams of site- 
centered dark breather solutions are shown in Fig. 12, which suggests that the site-centered 
dark breather solutions with frequency close to the linear frequency c o are long-lived and have 
marginal oscillatory instability, below the pertinent critical point. However, the oscillatory 
instability becomes more and more significant as oj b increases, leading to the breakdown of 
the dark breather structure of the solution. In fact, beyond the critical point, the breakup 
of the site-centered breather appears to be accompanied in Fig. [12] by a dramatic evolution, 
whereby the configuration is completely destroyed and a form of lattice dynamical turbulence 
ensues. This phenomenon is reminiscent of traveling wave instabilities observed in 31 for 


Hertzian chains and may be worth further study, which, however, is outside the scope of the 
present manuscript. 

In contrast to the site-centered solutions, the bond-centered ones exhibit only real insta¬ 
bility at the early stage of the continuation when the breather frequency c o b is greater than 
but close to uj. At those frequencies, the magnitude of the Floquet multipliers corresponding 
to the real instability of bond-centered breathers is larger than the moduli of the multipliers 
describing the oscillatory instability of the site-centered ones, resulting in not only shorter 
lifetime of the bond-centered solutions, but also setting the dark breather state in motion. 
A representative space-time evolution diagram for bond-centered dark breather solution of 


frequency oj b = 2.05 is shown in the left panel of Fig. 13 In the right panel of Fig. 13 


the manifestation of real instability of the same solution is shown, where the perturbation 
of the dark breather solution along the direction associated with the unstable mode that 
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Figure 11: Left plot: moduli of the Floquet multipliers of the site-centered breathers for frequency Ub < 
2.08. Right plot: the relative error Eb(t) versus frequency at t = 1200. Inset shows the relative error for 
frequencies less than ujb < 2.077. Here n = 1 and p = 1/3. 



Figure 12: Left panel: contour plot of the time evolution of the site-centered solution for wj, = 2.05. 
The color bar corresponds to the magnitude of the strain x n (top) and y n (bottom). Right panel: same 
computation as in left panel but for cob = 2.09. Here k = 1 and p = 1/3. 
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corresponds to a real Floquet multiplier is used as the initial condition for the integration. 
One can see that the instability results in a dark breather moving with constant velocity 
after some initial transient time in the left panel, while the pertinent motion is initiated 
essentially immediately by the perturbation induced in the right panel. However, as LOb is 



Figure 13: Left panel: contour plot of the time evolution of the bond-centered solution for ojb = 2.05. 
The color bar corresponds to the magnitude of the strain x n (top) and y n (bottom). Right panel: same 
simulation as in the left panel but with the perturbed dark breather as the initial condition. Here u = 1 and 
P = 1/3- 


increased the same phenomenology (dismantling of the breather and chaotic evolution) is 
also taking place for the bond-centered breathers, as shown in Fig. [Td| 



Figure 14: Left panel: contour plot of the time evolution of the bond-centered solution for u>b = 2.08. 
The color bar corresponds to the magnitude of the strain x n (top) and y n (bottom). Right panel: same 
simulation as in the left panel but at frequency u>b = 2.10. Here k = 1 and p = 1/3. 


The large arc seen in the middle of the Floquet multiplier diagram of the bond-centered 
breather solutions (the top left plot in Fig. 10) corresponds to the period-doubling bifurcation. 
As the frequency approaches u^ ~ 2.063, two complex conjugate eigenvalues collide on the 
real axis at —1. Two real eigenvalues then form and move in opposite directions as uJb 
increases. After the difference between the real eigenvalues reaches a maximum value, they 

























































start moving toward each other and collide at oo b ~ 2.096. Breathers with double the period 
(half the frequency) of the ones on the main branch exist between these two frequencies. 

To explore the numerically exact period-doubling orbits, we constructed an initial seed 
by slightly perturbing the dark breather solution at the bifurcation point along the direction 
of eigenvector associated with the eigenvalue —1. As shown in Fig. 15, the eigenvector is 
spatially localized at the middle of the chain. As a consequence, the initial seed only differs 
from the previous dark breather solution in the middle part of the chain. Sample profiles of 



Figure 15: Left panel: the eigenvector associated with the eigenvalue —1. Right panel: numerically exact 
dark breather solution of frequency oj b = 1.032 (connected squares) and the initial seed (connected stars) 
obtained by perturbing numerically exact dark breather solution of frequency u] b = 2.063. 


period-doubling dark breather solution at tu b = 1.032 are shown in the top of Fig. 16 To check 


that these solutions differ from the main branch, we integrated the solution for both the full 
period T b = 2 tt/ oj h and its half T b /2 and verified that the period of the obtained new solution 
is doubled compared with the previously obtained dark breathers. The continuation method 
is used to obtain all the period-doubling dark breather solutions for different frequencies. 
Note that the continuation stops at co b ~ 1.051 and also cannot proceed for co b below 1.032, 
which agrees with the (doubled) frequency range of the large arc in the top left plot of 
Fig. 10 All of these solutions exhibit both real and oscillatory instabilities (see the bottom 


plots of Fig. |16J). 

We now investigate the effect of mass ratio p on the system by repeating the same exper¬ 
iment for relatively large and small values p. Recall that the linear frequency uj = \Jk + n/p 
decreases and approaches A J~k as p increases. In what follows, we start the continuation at fre¬ 


quency u b = (u+0.01 to make sure that the amplitude of the initial seed (89) is small. We first 
test the case when p = 3 with linear frequency given by uj — 1.1547. Sample profiles of bond- 
centered and site-centered dark breather solutions at the frequency co b = lo + 0.1 = 1.2547 
and their space-time evolution diagrams are shown in Fig. [TT} We observe the bond-centered 
solution starts to move in form of a traveling dark breather after the integration for a suffi¬ 
ciently long time, while the site-centered solution persists for a longer time in the simulation 
and hence can be considered to be long-lived. The dynamic behavior of both solutions is 
consistent with their numerically computed Floquet spectrum shown in the right of Fig. [I8j 


Moreover, the diagrams of Floquet multipliers’ moduli in the right panel of Fig. 18 suggest 
that the bond-centered solutions exhibit only real instability for a wide range of frequencies 
[uj + 0.001, uj + 0.2], while the site-centered dark breathers have just marginal oscillatory 
instability. 

Next, we repeat the simulation with p = 10 and c o = 1.0488. Note that the amplitude of 
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Figure 16: Top panel: comparison of the numerically exact double-period solution at u>b = 1.032 (connected 
squares) and the result of its integration (connected stars) at t, = 71/2 and t = T b . Bottom left: moduli 
of Floquet multipliers versus frequency uib for the period-doubling solutions. The Floquet multipliers for 
LOb = 1-032 in the complex plane are shown in the right plot. Here k = 1 and p = 1/3. 




Figure 17: Top panel: sample profiles (circles) of site-centered (left) and bond-centered (right) dark 
breathers at the frequency ojb = 1.2547. Stars connected by dashed lines represent strain profiles after the 
integration over 299Tb « 1500. Note that the site-centered solution has relative error F?{,(1500) = 9.21 x 10 -5 . 
Bottom panel: space-time evolution diagrams for site-centered (left) and bond-centered (right) solutions. 
Here k = 1 and p = 3. 
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Figure 18: Left panel: moduli of Floquet multipliers versus frequency Wf, for the bond-centered (top) 
and site-centered (bottom) types. Right panel: Floquet multipliers of dark breather solutions of frequency 
u>b = 1.2547 in the complex plane. Here k = 1 and p = 3. 


dark breather solution at frequency uij, = cu + 0.01 is close to 1.6 x 10 4 . As shown in Fig. 19 


the pattern of Floquet multipliers moduli is very similar to the p = 3 case for breather 
frequencies close to u. However, as u>b becomes larger, we observed significant oscillatory 
instability of both the bond-centered and site-centered solutions. Note that the distribution 
of Floquet multipliers at large mass ratios (for example, p = 3, 10) is completely different 
than that at smaller ones (such as p = 1/3), given the same magnitude of frequency difference 

OJb — OJ. 

We now fix the breather frequency cuf, and perform the continuation in the mass ratio p. 


The results of numerical continuation are shown in Fig. 20 At a given breather frequency, the 
real instability is only exhibited by solutions of the bond-centered type, and its significance 
is gradually increasing as mass ratio becomes larger. In contrast, the site-centered solutions 
have marginal oscillatory instability and persist for a long time. Moreover, the lifetime of 
those solutions decreases as the mass ratio increases. At a larger frequency like u>b = 2.05 and 
small mass ratio, the emergence of many unstable quartets suggests that both bond-centered 
and site-centered solutions share strong modulational instabilty of the background, which 
leads to a chaotic evolution of both solutions after a short time of integration. 


8 Concluding remarks 

In this work, we studied nonlinear waves in a resonant granular material modeled by a 
Hertzian chain of identical particles with a secondary mass attached to each bead in the 
chain by a linear spring. Following the approach developed in |[l5) for a limiting case of 
the present model, we derived generalized modulation equations of DpS type. We showed 
that for suitable initial data and large enough mass ratio, these equations reduce to the DpS 
equation derived in [15) and rigorously justified the equation in this limit on the long time 
scales. We then used the DpS equations to investigate the time-periodic traveling wave of 
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Figure 19: Left panel: moduli of Floquet multipliers versus frequency Wb for the bond-centered (top) 
and site-centered (bottom) types. Right panel: Floquet spectrum of dark breather solutions of frequency 
uib = 1-20 in the complex plane. The inset represents the zoom-in of moduli of Floquet multipliers for the 
site-centered solution at the frequencies Ub G [w + 0.001, w + 0.07]. Here k = 1 and p = 10. 



Figure 20: Moduli of Floquet multipliers versus mass ratio p for the bond-centered (top) and site-centered 
(bottom) type. The tested frequencies of the dark breather are ujb = 1-50,1.75 and 2.05. Here k = 1. 
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the system at finite mass ratio. We showed numerically that these equations can successfully 
capture the dynamics of small-amplitude periodic traveling waves. 

Turning our attention to the breather-type solutions, we proved non-existence of nontriv¬ 
ial bright breathers at finite mass ratio. However, we also showed that at sufficiently large 
mass ratio and suitable initial data, the problem has long-lived bright breather solutions. 

The generalized DpS equations were also used to construct well-prepared initial condi¬ 
tions for the numerical computation of dark breather solutions. A continuation procedure 
based on a Newton-type fixed point method and initiated by the approximate dark breather 
solutions obtained from the DpS equations was utilized to compute numerically exact dark 
breathers for a wide range of frequencies and at different mass ratios. The stability and 
the bifurcation structure of the numerically exact dark breathers of both bond-centered and 
site-centered types were examined. Our numerical results strongly suggest that the bond- 
centered solutions exhibit real instability that may give rise to steady propagation of a dark 
breather after large enough time. In addition, period-doubling bifurcations of these solutions 
were identified at small mass ratios. The site-centered solutions, in contrast to the bond- 
centered ones, appeared to exhibit only oscillatory instability, which is much weaker than 
the real instability of the bond-centered breathers for a range of breather frequencies that 
are close enough to the natural frequency of the system, i.e., the frequency of out-of-phase 
motion within each unit cell of the chain involving the particle and the secondary mass. As 
a consequence, these low-frequency site-centered solutions persisted for a long time in the 
numerical simulation, and thus the effect of oscillatory instability is quite weak. However, 
we also provided case examples of their (long-time) instabilities that led to their complete 
destruction and ensuing apparently chaotic dynamics within the lattice. We showed that 
the distribution of Floquet multipliers and hence stability of the dark breather solutions are 
significantly affected by the mass ratio and breather frequency. 

A challenge left for the future work is to rigorously prove the existence of small-amplitude 
exact periodic traveling wave and dark breather solutions of system (J2]) using the approxi¬ 
mate solutions obtained from the generalized DpS equations. Another intriguing aspect to 
further consider involves the mobility of the dark breathers, and its association with the 
dynamical instability of the states, as well as possibly with the famous Peierls-Nabarro bar¬ 
rier associated with the energy difference between bond- and site-centered solutions, i.e., the 
energy barrier that needs to be “overcome” in order to have mobility of the dark breathers. 
Equally important and relevant would be an effort to analytically understand the modula- 
tional stability properties of the lattice, perhaps at the DpS level and compare them with 
corresponding systematic numerical computations. On the experimental side, it will be in¬ 
teresting to investigate whether we can generate dark breathers by exciting the both ends of 
a finite chain in a way similar to [7|. Additionally, exciting small amplitude traveling waves 
through boundary excitations, e.g. in the woodpile chain of 22 and observing experimen¬ 


tally their evolution through lased Doppler vibrometry would also be particularly relevant. 
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